#include "exportToVTK.h"

void exportSolutionToVTK(pMesh theMesh, void *pData1, void *pData2, void *pData3, string filename){
	// open file
	ofstream fid;
	fid.open(filename.c_str());

	// check if it's OK, otherwise terminate program
	char msg[512]; sprintf(msg,"File '%s' could not be opened or it does not exist.\n"
			"\tCheck if directory was typed correctly.\n",filename.c_str());

	// print data to file
	fid << "# vtk DataFile Version 2.0\n";
	fid << "Two phases flow simulation\n";
	fid << "ASCII\nDATASET UNSTRUCTURED_GRID\n";
	fid << "POINTS " << M_numVertices(theMesh) << " float\n";

	int dim = theMesh->getDim();
	int numElements = getNumElements(theMesh,dim);
	printVerticesCoordenates(fid,theMesh);
	printElementConnectivities(fid,theMesh,dim,numElements);
	printCellTypeList(fid,dim,numElements);
	fid.close();
}

// print vertex coordenates and transfer to each one all computed values
// to be printed later
// = = = = = = = = = = = = = = = =  = = = = = = = = = = = = = = = = = =
void printVerticesCoordenates(ofstream &fid, pMesh theMesh){
	pEntity e;
	int count = 0;
	VIter vit = M_vertexIter(theMesh);
	while( (e = VIter_next(vit)) ) {
		double coord[3] = {.0, .0, .0};
		V_coord(e,coord);
		for (int i=0; i<3; i++) fid << coord[i] << " ";
		fid << endl;
		EN_attachDataInt(e,MD_lookupMeshDataId("mLN"),count++);
	}
	VIter_delete(vit);
}

// print elements connectivities
void printElementConnectivities(ofstream &fid, pMesh theMesh, int dim, int numElements){
	fid << "\nCELLS " << numElements << " " << (dim+2)*numElements << endl;
	pEntity elem;
	if (dim==2){
		FIter fit = M_faceIter(theMesh);
		while(elem = FIter_next(fit))
			if ( !theMesh->getRefinementDepth(elem) )
				printElementConnectivities(fid,elem,dim);
		FIter_delete(fit);
	}
	else{
		RIter rit = M_regionIter(theMesh);
		while (elem = RIter_next(rit))
			if ( !theMesh->getRefinementDepth(elem) )
				printElementConnectivities(fid,elem,dim);
		RIter_delete(rit);
	}
}

void printElementConnectivities(ofstream &fid, pEntity elem, int dim){
	int mappedLNodes;
	fid << dim+1 << " ";
	for(int i=0; i<dim+1; i++){
		EN_getDataInt((pVertex)elem->get(0,i),MD_lookupMeshDataId("mLN"),&mappedLNodes);
		fid << mappedLNodes << " ";
	}
	fid << endl;
}

void printCellTypeList(ofstream &fid, int dim, int numElements){
	fid << "\nCELL_TYPES " << numElements << endl;
	int type = (dim==2)?5:10;
	for(int i=0; i<numElements; i++)
		fid << type << endl;
}

int getNumElements(pMesh theMesh, int dim){
	int numElements = 0;
	pEntity elem;
	if (dim==2){
		FIter fit = M_faceIter(theMesh);
		while(elem = FIter_next(fit))
			if ( !theMesh->getRefinementDepth(elem) )
				numElements++;
		FIter_delete(fit);
	}
	else{
		RIter rit = M_regionIter(theMesh);
		while (elem = RIter_next(rit))
			if ( !theMesh->getRefinementDepth(elem) )
				numElements++;
	}
	return numElements;
}


